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ABSTRACT 


System identification is utilised in the aerospace community for development of simulation 
models for robust control law design. These models are often described as linear, time-invariant 
processes and assumed to be uniform throughout the flight envelope. Nevertheless, it is well known 
that the underlying process is inherently nonlinear. Over the past several decades the controls and 
biomedical communities have made great advances in developing tools for the identification of 
nonlinear systems. 

In this report, we show the application of one such nonlinear system identification technique, 
structure detection, for the analysis of Quiet Spike™ (Gulfstream Aerospace Corporation, Savannah, 
Georgia) aeroservoelastic flight-test data. Structure detection is concerned with the selection of a 
subset of candidate terms that best describe the observed output. Structure computation as a tool 
for black-box modeling may be of critical importance for the development of robust, parsimonious 
models for the flight-test community. 

The objectives of this study are to demonstrate via analysis of Quiet Spike™ aeroservoelastic 
flight-test data for several flight conditions that: linear models are inefficient for modelling 
aeroservoelastic data, nonlinear identification provides a parsimonious model description whilst 
providing a high percent fit for cross-validated data and the model structure and parameters vary 
as the flight condition is altered. 


NOMENCLATURE 


AIC 

LTI 

MDL 

NARMAX 

NL 


Akaike’s information criterion 
linear time-invariant 
minimum description length 

Nonlinear AutoRegressive, Moving Average exogenous 
nonlinear 


INTRODUCTION 

System identification, or black-box modeling, is a critical step in aircraft development, 
analysis and validation for flightworthiness. The development and testing of aircraft typically takes 
many years and requires considerable expenditure of limited resources. One reason for lengthy 
development time and cost is assuming the underlying system is linear and invariant throughout 
the flight envelope. This assumption is related to having an inadequate knowledge of an appropriate 
model type or structure to use for parameter estimation. Selection of an insufficient model structure 
may lead to difficulties in parameter estimation, giving estimates with significant biases and/or 
large variances (ref. 1). This often complicates control synthesis or renders it infeasible. The power 
of using structure detection techniques as a tool for model development (i.e., black-box modeling) 
is that it can provide a parsimonious system description which can describe complex aeroelastic 


behaviour over a large operating range. Consequently, this provides models that can be more robust 
and, therefore, reduce development time. 

Moreover, when studying aeroelastic systems it may not be practical to assume that the exact 
model structure is well known a priori. In aerospace systems analysis one of the main objectives is 
not only to estimate system parameters but also to gain insight into the structure of the underlying 
system. Therefore, structure computation is of significant relevance and importance to modelling 
and design of aircraft and aerospace vehicles. Structure computation may indicate deficiencies in an 
analytical model and could lead to improved modelling strategies and also provide a parsimonious, 
black-box, system description for control synthesis (ref. 2). 

For linear systems modelling a commonly used approach for determining model structure is 
the minimum description length (MDL) proposed by Rissanen (ref. 3). The MDL was specifically 
developed to overcome some of the inconsistencies of the Akaike’s information criterion (AIC) 
(ref. 4), i.e., its variance does not tend to zero for larger sample sizes, N. 

Recently, the bootstrap method has been shown to be a useful tool for structure detection 
of nonlinear models (refs. 5-7). The bootstrap is a numerical method for estimating parameter 
statistics which requires few assumptions (ref. 8). The conditions needed to apply bootstrap to least- 
squares estimation are quite mild; namely, that the errors be independent, identically distributed, 
and have zero-mean. 

Identification of flight dynamics is not well understood in the subsonic, transonic and 
supersonic regimes. Therefore, in this report, we investigate whether (1) a linear or nonlinear 
model best represents the observed data and (2) the system is invariant during envelope expansion 
(varying Mach number). The data analysed in this study is from the F-15B (McDonnell Douglas, 
now the Boeing Company, Chicago, Illinois) Quiet Spike™ (Gulfstream Aerospace Corporation, 
Savannah, Georgia) flight-test program which was a collaborative effort between Gulfstream 
Aerospace Corporation and the NASA Dryden Flight Research Center (refs. 9-11). The results 
show that (1) linear models are inefficient for modelling aeroservoelastic data, (2) nonlinear 
identification provides a parsimonious model description whilst providing a high percent fit for 
cross- validated data and (3) the model structure and parameters vary as the flight condition is 
altered. 


NARMAX MODEL FORM 

The dynamic behavior of many nonlinear systems can be represented as a discrete polynomial 
which expands the present output value in terms of present and past values of the input signal 
and past values of the output signal (refs. 12- 14). A system modelled in this form is popularly 
known as a NARMAX (Nonlinear AutoRegressive, Moving Average exogenous) model and is 
linear- in-the -parameters. 
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Recently, Kukreja et al. (ref. 15) showed that NARMAX identification is well suited to 
describing aeroelastic phenomena. The NARMAX structure is a general parametric form for 


modelling nonlinear systems (ref. 12). This structure describes both the stochastic and deterministic 
components of nonlinear systems. Many nonlinear systems are a special case of the general 
NARMAX structure (ref. 16). In this report, we focus on a special class of NARMAX models: 
nonlinear polynomial models. The polynomial NARMAX structure models the input-output 
relationship as a nonlinear difference equation of the form shown in eq. (1): 


/denotes a nonlinear mapping, 1 is the order of the nonlinearity, u is the controlled or exogenous 
input, z is the measured output, and e is the uncontrolled input or innovation. This nonlinear mapping 
may include a variety of nonlinear terms, such as terms raised to an integer power, products of 
present and past inputs, past outputs, or cross-terms. In general, the nonlinear mapping ,/ can be 
described by a wide variety of nonlinear functions such as sigmoids or splines (refs. 17, 16). This 
system description encompasses many forms of linear and nonlinear difference equations that are 
linear- in-the -parameters . 

Identifying a NARMAX model requires two things: (1) structure detection and (2) parameter 
estimation. Structure detection can be divided into: (la) model order selection and (lb) selecting 
which parameters to include in the model. We consider model order selection as part of structure 
detection since, theoretically, there are an infinite number of candidate terms that could be considered 
initially. Establishing the model order, then, limits the choice of terms to be considered. For the 
NARMAX model, the system order is defined to be an ordered tuple as shown in eq. (2): 


where n u is the maximum lag on the input, n z the maximum lag on the output, n e the maximum 
lag on the error and l is the maximum nonlinearity order. Note that for non-polynomial NARMAX 
models, l may be simply replaced by a nonlinear mapping of some specified class. In this report, 
we assume that the system order is known. 


The structure detection problem is that of selecting the subset of candidate terms that 
best describes the output. Therefore, the parametrisation of a system is still further reduced by 
determining which of the components are required. The maximum number of terms in a NARMAX 
model with n z , n u and n e dynamic terms and Zth order nonlinearity is shown in eq. (3): 



0 = [n u n z n e l] 


( 2 ) 


STRUCTURE DETECTION 



( 3 ) 
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As a result, the number of candidate terms becomes very large for even moderately complex 
models making structure detection difficult. We define the maximum number of terms, p, as 
the number of candidate terms to be initially considered for identification. Due to the excessive 
parameterisation (the curse of dimensionality), the structure detection problem often leads to 
computationally intractable combinatorial optimisation problems. 

TIME SERIES 

The data considered in this report is time series since the input signal, u(n), is assumed to be 
zero or constant. Time series analysis is used when inputs are not available to the experimenter, or 
where it is unclear which signals are inputs and which are outputs (ref. 18). Models arising from 
time series data can have many forms (ref. 1). However, in our treatment of the data the ARMA 
and NARMA model class are of practical significance. 

This special case of the general NARMAX model [eq. (1)] can be written as shown in 
eq. (4): 


z(n) = f l \z(n-l),...,z(n-n z )...,e(n-n e )l + e(n) 


( 4 ) 


where we redefine the model order for this model set as shown in eq. (5): 

O = 


( 5 ) 


The maximum number of candidate terms in a model [eq. (4)] with n z and 
and /th order nonlinearity is shown in eq. (6): 


n e dynamic terms 


P = 


l 

n 

i = 1 


n z + n e + 1 


(6) 


Note that ARMA models can be estimated using the Yule-Walker equations or the instrumental 
variable (IV) estimator to avoid estimating the MApart (ref. 1). This is the approach taken in this 
report. For a NARMA model the NMA part must be modelled. For nonlinear systems, output 
additive noise can produce multiplicative terms between input, output and itself. To compute 
unbiased parameter estimates a noise model (i.e. NMA) needs to be estimated (ref. 19). 

STRUCTURE DETECTION METHODS 

With the model types defined for the flight-test data available for analysis, we describe two 
approaches applicable to these model classes. The first is appropriate for AR models whilst the 
second is appropriate for NARMA models. 
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Minimum Description Length 


A commonly used technique in linear system identification to determine model structure is 
MDL (ref. 3). The MDL was specifically developed to overcome some of the inconsistencies of 
the Akaike’s information criterion (AIC) (ref. 4), i.e. its variance does not tend to zero for larger 
sample sizes, N. 

The number of parameters necessary to reproduce an observed sequence { Z\ , . . . , Z /V } of a 
time series depends on the model and parameters assumed to have generated the data (ref. 3). The 
MDL technique finds the model which minimises the description length and thereby computes an 
estimate of model order (ref. 3). 

Binary prefix codes are used to encode data strings. These data strings can be made up of 
symbols, parameters, numbers, etc. It is known that the average length of a code word is bounded 
by Shannon’s theorem (ref. 3). Therefore, it is possible to write eq. (7) (ref. 3): 

(x) L (x) >-£/?(*) log p(x) 

X X 

where L(x) is the length of the code word (i.e. length of parameter vector 9) and p(x) is the 
probability of x. It is also possible to write eq. (8): 

L{z\x,e) = -\ogp{z\x,d) (8) 


where L(z\x, 9 ) is known as the log-likelihood function (to be maximised). Let 9 denote the value of 
the parameter which maximises the likelihood and thus minimises the parameter vector length (i.e. 
code word length) L{z\x, 9 ). Since 6 can only be encoded up to a certain precision, the code word 
length, L(z\x, 9), becomes longer than the desired minimum L(z\x, 9 ), given noise considerations. 
Let the precision be S = 2~ q where q is the number of bits used for encoding the parameter. It is 
possible to save on the code word length if q is small. However, the result is a loss in precision. 
The optimal precision depends on the size of the observed data via - log 8 = 0.5 log N , and hence the 
total code word length for k parameters is given by the MDL, eq. (9): 


MDL(k) 


log [maximised likelihood] + ^ k log N 


which, for an AR(« Z ) model gives eq. (10): 


(9) 


MDL(«,) = log[maximised likelihood] + —log A. 


( 10 ) 


5 


Bootstrap 


Recently, the bootstrap has been shown to be a useful tool for structure detection of nonlinear 
models (ref. 6). The bootstrap is a numerical method for estimating parameter statistics which 
requires few assumptions (ref. 8). The conditions needed to apply bootstrap to least-squares 
estimation are quite mild; namely, that the errors be independent, identically distributed, and have 
zero-mean. 


The bootstrap is a technique to randomly reassign observations which enables re-estimates of 
parameters to be computed. This randomisation and computation of parameters is done numerous 
times and treated as repeated experiments. In essence, the bootstrap simulates a Monte Carlo 
analysis. For structure computation, the bootstrap method is used to detect spurious parameters, 
those parameters whose estimated values cannot be distinguished from zero. 


Application of an appropriate £2 estimator to measured data gives the model response Z, 
residuals £ and parameter estimate 6 .The bootstrap assumes that these residuals, e = [ej, £ 2 , ] , 
arise from an unknown distribution, D. By randomly resampling these residuals, with replacement, 

it is possible to generate a resampled version of the prediction errors, e* = ej ,£ 2 ,---£n , whose 

distribution estimates D. The resampling procedure for each i* involves randomly selecting from 
e with an equal probability associated with each of the N elements. For example, a possible 

resampled version of the errors for N = 5 is £* = [£ 4 , ^,£ 4 , £ 2 , £ 3 ]. The star notation indicates that 
£* is not the calculated error £ but rather a resampled version of it. These resampled errors are 


added to the model response to generate a bootstrap replication of the original data, shown in 
eq. (11): 


z ='v ze e+£ . 


(ll) 


A new bootstrapped parameter estimate 6 is obtained from this bootstrapped data Z ' . 
This procedure is repeated B times to provide a set of parameter estimates from the B bootstrap 
replications, as shown in eq. (12): 


0 * = 


9\,...,9 B 


( 12 ) 


Parameter statistics can then be easily computed from 0 by forming percentile intervals at 
a chosen level of significance, a . 

Structure detection can provide useful process insights that can be used in subsequent 
development or refinement of physical models. Therefore, in the sequel, we investigate the 
applicability of MDL and the bootstrap to experimental aircraft data. Specifically, MDL and 
bootstrap methods are used as structure detection tools to assess whether the (1) underlying data 
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is best described by a linear time-invariant (LTI) or nonlinear model and (2) model structure is 
invariant during envelope expansion. 

EXPERIMENTAL F-15B QUIET SPIKE™ DATA 

The MDL and bootstrap technique were applied to experimental flight-test data from the 
F-15B Quiet Spike™ project by Gulfstream Aerospace Corporation and the NASA Dryden Flight 
Research Center (refs. 9-11). The data analysed for this study used structural accelerometer 
response output of the Quiet Spike™ boom tip when fully extended. 

Procedures 

Flight data was gathered during subsonic, transonic and supersonic flutter clearance of the 
F-15B Quiet Spike™ (fig. 1). This report considers accelerometer data measured during pilot 
induced pitch-raps at Mach 0.85, 0.95 and 1.40 at an altitude of 12,192 m (40,000 ft). The output 
was taken as the response of an accelerometer mounted near the Quiet Spike™ boom tip (fig. 2). 
Data was sampled at 400 Flz. For analysis, the recorded flight-test data was decimated by a factor 
of 2, resulting in a final sampling rate of 200 Hz. 



ED06-0 184-33 

Figure 1. Flight test article in extended configuration. 
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Primary 



Figure 2. Quiet Spike™ sensor location. 


For identification of a linear model an arbitrarily large AR model of fiftieth order ( n z = 50) 
was posed for estimation and the MDL technique used to determine the optimal linear model. For 
identification of a nonlinear model a model order with fourth-order output and error dynamics 
and second-order nonlinearity, O = [4, 4, 2], was used. A model with fourth-order dynamics was 
selected because it has been observed that aeroservoelastic structures are well defined by a fourth- 
order LTI system (ref. 20). The nonlinearity order was chosen as second-order because empirical 
results showed models of higher nonlinear order were not efficient to describe the data sets 
available for analysis. This gave a full model description with 45 candidate terms. The nonlinear 
model was identified applying the bootstrap approach. For the bootstrap method, B = 100 bootstrap 
replications were generated to assess the distribution of each parameter. For all three techniques, 
each parameter was tested for significance at the 95 percent confidence level. 

For both linear and nonlinear identification, table 1 shows the number of data points available 
for each flight condition. The estimation data sets were from accelerometer response measurements 
of the primary sensor on the boom tip and the cross-validation data sets were of the backup sensor 
at the same location. 


Table 1. Data points available at each flight condition. 


Altitude 12,192 m 

40,000 ft) 

Mach number 

0.85 

0.95 

1.40 

Estimation data 

572 

400 

504 

Validation data 

572 

400 

504 
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Results 


The results of identifying the F- 1 5B Quiet Spike™ data are presented. Figures 3(a)-3(c) show 
the output data sets used for this analysis. The data represents structural accelerometer response 
(primary sensor boom tip) used to compute the system structure. 



(a) Mach 0.85. 

Figure 3. Estimation data: recorded structural accelerometer response to stick-raps. 


9 



(b) Mach 0.95. 



(c) Mach 1.40. 
Figure 3. Concluded. 
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Equations (13)— (15) depict the model structure computed by the bootstrap method. 


Mach 0.85, altitude 40,000 ft (12,192 m) 

z(n ) = 6\z(n - \) + Q2z{n-2) + OT,z{n-'i) + 0^ {n-4) + G$£(n-\) 

+d(,e(n -2)+Gj£(n- 3) + 0gz(« — 4)e(n — 4) + Gg£~ (n- 4) + e[n) 

Mach 0.95, altitude 40,000 ft (12,192 m) 
z(n) = Y\z(n - 1)+ 2) + 73 z(/z - 3)+ 74 z(n- 3 )z(n- 4) 

+f 5 Z 2 (n ~ 4) + y 6 e(n - l) + f 7 e(n - 2) + y 8 £ (n - 3) 

+Ygz(n — 3)e(/i -4) + Y\qz(h - 4)s(«- 3) + 7 ^e(n- 3)e(n- 4) 
+ 7 l 2 z(« - 4 )e(n - 4) + y 13 £ 2 (n - 4 ) + e(n) 

Mach 1.40, altitude 40,000 ft (12,192 m) 

z(n) = P\z(n - 1) + P 2 z ( n - 2) + Py.(n - 3) + ji^z{n - 1 )z(n - 4) 

+/I 5 £(n - 1) + P(,£(n-2) + e(«- 3) + P 8 z(n- l)£(n- 4) 

+Pg z(n — 4)e(n — l) + Piq£(h - 1 )s(« - 4) + e(n). 


( 13 ) 


(14) 


(15) 


Equations 13 through 15 represent the computed model structure for flight conditions Mach 
0.85, 0.95 and 1.40, respectively. The computed model structures are represented by a combination 
of linear and nonlinear, lagged input-output terms and contain 9, 13 and 10 terms for Mach 0.85, 
0.95 and 1.40, respectively. Hence, the bootstrap technique successfully produced a parsimonious 
model description from the full set of 45 candidate terms. 

For AR (linear) model identification using MDL to compute structure, the estimated models 
were of order n z = 42, 44, and 46 for Mach 0.85, 0.95 and 1.40, respectively. These models are not 
shown since they are simply a dynamic expansion of the output up to the order stated. However, 
for cross-validation data, we show the model fit of these linear models compared to the cross- 
validation fit obtained with the NARMA models [see figs. 4(a)-4(c)]. 
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(a) Mach 0.85. 


Cross-validated structural accelerometer response 



(b) Mach 0.95. 

Figure 4. Cross-validation data: predicted linear and nonlinear model accelerometer response of 
z-tip longitudinal sensor superimposed on measured accelerometer output. 
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(c) Mach 1.40. 
Figure 4. Concluded. 


Figures 4(a)-4(c) show the predicted output for the cross-validation data sets for the identified 
structures [(a): eq. (13) and AR(«, = 42); (b): eq. (14) and AR( n z = 44); and (c): eq. (15) and 
AR(n„ = 46)]. Each panel displays the full time history of the predicted output of the linear and 
nonlinear models superimposed on top of the measured output. For Mach 0.85 [fig. 4(a)] the 
predicted output of the linear and nonlinear models account for over 91 percent and 95 percent 
of the measured outputs variance, respectively. For Mach 0.95 [fig. 4(b)] the predicted output 
of the linear and nonlinear models account for over 92 percent and 97 percent of the measured 
outputs variance, respectively. For Mach 1.40 [fig. 4(c)] the predicted output of the linear and 
nonlinear models account for over 91 percent and 96 percent of the measured outputs variance, 
respectively. 

The results demonstrate that although the AR models contain many more terms to explain 
the underlying process they still offer a lower percent fit compared to the nonlinear model at the 
cost of model complexity (higher order) which often leads to more complex control synthesis. 
The nonlinear models contain only a few terms and were capable of explaining a larger percent of 
the output variance. For these data sets the results show linear models are inefficient for accurate 
modelling of aeroservoelastic data. These results show a nonlinear identification approach offers 
a parsimonious system description whilst providing a high percent fit for cross-validated data. 
Moreover, the results illustrate the need to vary model structure for different flight conditions. 
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DISCUSSION 


Experimental results demonstrate that structure computation as a tool for black-box modeling 
may be useful for the analysis of dynamic aircraft data. The bootstrap successfully reduced the 
number of regressors posed to aircraft aeroelastic data yielding a parsimonious model structure for 
each data set. Additionally, these parsimonious structures were capable of predicting a large portion 
of the cross-validation data, collected on a backup sensor at a similar location. However, for linear 
analysis, the MDL approach was not able to reduce the model order (structure) as well and yielded 
a more complex system description. Whilst these linear models have higher complexity (degrees 
of freedom), they provided a model predictive capability that explained a smaller percent of the 
observed output variance. This find indicates a linear model may not be appropriate to describe 
aeroservoelastic data. A higher percent fit offered by the parsimonious nonlinear models suggests 
that the identified structures and parameters explain the data well. Using percent fit alone as an 
indicator of model goodness could lead to incorrect interpretations of model validity. Nevertheless, 
in many cases for nonlinear models this may be the only indicator that is readily available. 

In this work, the results show that whilst the linear dynamics remained invariant for all flight 
conditions available for analysis, the nonlinear dynamics changed as Mach number increased. For 
Mach 0.85 the model [eq. (13)] displayed a mildly nonlinear process which physically makes sense 
since the airflow is mainly subsonic. When the Mach number was increased to 0.95 the model 
[eq. (14)] demonstrated a richer nonlinear dynamic description which is likely due to embedded 
shock formations in the transonic regime. For Mach 1.40 the model [eq. (15)] displayed a mildly 
nonlinear process again which physically makes sense since in this regime the shocks become 
fixed. It is difficult to make definitive comments on the underlying physics responsible for this 
behaviour without extensive analysis of different flight conditions. The important points to note 
are, this study suggests (1) nonlinear models are appropriate to describe the dynamics behaviour 
of advanced aircraft and (2) models describing aircraft dynamics vary with flight condition. This 
suggests nonlinear modelling may afford a robust and parsimonious system description over a 
larger operating regime and models used for prediction (e.g. control) should not be invariant for all 
flight conditions. This may hold significant implications for aircraft development. 

For this study, only a polynomial mapping with fourth-order output and error lag was used 
as a basis function to explain the nonlinear behaviour of the F-15B Quiet Spike™ data. Clearly, 
different basis functions and a higher dynamic order (lag-order) should be investigated to determine 
if another basis can produce accurate model predictions with reduced complexity. Moreover, further 
studies are necessary to evaluate whether the model structure is invariant under different operating 
conditions, such as altitude, and model parameterisations. 

This study illustrates the usefulness of structure detection as an approach to compute a 
parsimonious model of a highly complex nonlinear process, as demonstrated with experimental 
data of aircraft aeroelastic dynamics. Moreover, analysis of flight-test data can provide useful 
process insights that can be used in subsequent development or refinement of physical models. 
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In particular, morphological models are based on assumptions (e.g. these effects are important 
and those are negligible) which may be incorrect (refs. 2 1 and 22). A structure computation approach 
to model identification may help uncover such surprises. 

CONCLUSIONS 

Results show that linear models are inefficient for modelling aeroservoelastic data and 
nonlinear identification provides a parsimonious model description whilst providing a high percent 
fit for cross-validated data. Moreover, the results demonstrate that model structure and parameters 
vary as the flight condition varies. These results may have practical significance in the analysis of 
aircraft dynamics during envelope expansion and could lead to more efficient control strategies. 
In addition, this technique could allow greater insight into the functionality of various systems 
dynamics, by providing a quantitative model which is easily interpretable. 
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